
/**
 * @file Error.h
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 
 * @version 0.1
 * @date 2020-12-02
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _ERROR_H_
#define _ERROR_H_
#include <Eigen/Dense>
#include <functional>
#include "Mesh.h"
#include "Element.h"
#include <iostream>
#define TEMPLATE template <unsigned int DIM>

TEMPLATE
class Error
{
    typedef Eigen::VectorXd VectorXd;
    typedef std::function<double(double*)> Function;
    
    public:
        Error(){std::cout << "construct!" << std::endl;};
        Error(Mesh<DIM> * mesh,Element<DIM> * element,Function real_solution_fun,VectorXd analytical_numerical_soluton_vec)
        {
            _mesh = mesh;
            _element = element;
            _real_solution_fun = real_solution_fun;
            _analytical_numerical_soluton_vec = analytical_numerical_soluton_vec;
        }
        ~Error(){};
        void setpartialderivativefunc(Function u_x,Function u_y){_real_solution_fun_x = u_x;_real_solution_fun_y = u_y;};
        double L2Norm();
        
        double LpNorm(double p);

        double H1Norm();
    private:
        Mesh<DIM> *_mesh;
        Element<DIM> *_element;
        Function _real_solution_fun;
        Function _real_solution_fun_x;
        Function _real_solution_fun_y;
        VectorXd _analytical_numerical_soluton_vec;
};

TEMPLATE   
double Error<DIM>::L2Norm()
{
    int n_Dofs = _element->n_Dofs();
    int n_GaussPnt = _element->n_GaussPnt();
    double l2error = 0.0;
    for(int k = 0;k < _mesh->n_element();k++)
    {
        double l2errorink = 0.0;
        std::vector<int> element_index = _mesh -> NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
	        temp[i] = _mesh ->DofsofIndex(element_index[i]);
        _element->SetDofsList(temp);
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j = 0.0;
            double u_h_j = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[DIM] = {x,y};
            u_j = _real_solution_fun(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j += _analytical_numerical_soluton_vec(element_index[i - 1]) * _element ->phi(xi,eta,i);
            l2errorink += (u_j - u_h_j) * (u_j - u_h_j) * _element ->GaussionWeight(j) * _element ->det_Jacobi(xi,eta);
        }
        l2error += l2errorink;
    }
    return sqrt(l2error);
}

TEMPLATE   
double Error<DIM>::LpNorm(double p)
{
    int n_Dofs = _element->n_Dofs();
    int n_GaussPnt = _element->n_GaussPnt();
    double l2error = 0.0;
    for(int k = 0;k < _mesh->n_element();k++)
    {
        double l2errorink = 0.0;
        std::vector<int> element_index = _mesh -> NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
	        temp[i] = _mesh ->DofsofIndex(element_index[i]);
        _element->SetDofsList(temp);
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j = 0.0;
            double u_h_j = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[DIM] = {x,y};
            u_j = _real_solution_fun(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j += _analytical_numerical_soluton_vec(element_index[i - 1]) * _element ->phi(xi,eta,i);
            l2errorink += pow((u_j - u_h_j),p) * _element ->GaussionWeight(j)*_element ->det_Jacobi(xi,eta);;
        }
        l2error += l2errorink;
    }
    return fabs(pow(l2error,1/p));
}

TEMPLATE   
double Error<DIM>::H1Norm()
{
    int n_Dofs = _element->n_Dofs();
    int n_GaussPnt = _element->n_GaussPnt();
    double h1error = 0.0;
    for(int k = 0;k < _mesh->n_element();k++)
    {
        double h1errorink = 0.0;
        std::vector<int> element_index = _mesh -> NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
	        temp[i] = _mesh ->DofsofIndex(element_index[i]);
        _element->SetDofsList(temp);
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j_x = 0.0;
            double u_h_j_x = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[DIM] = {x,y};
            u_j_x = _real_solution_fun_x(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j_x += _analytical_numerical_soluton_vec(element_index[i - 1]) * _element ->phi_x(xi,eta,i);
            h1errorink += (u_j_x - u_h_j_x) * (u_j_x - u_h_j_x) * _element ->GaussionWeight(j);
        }
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j_y = 0.0;
            double u_h_j_y = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[DIM] = {x,y};
            u_j_y = _real_solution_fun_y(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j_y += _analytical_numerical_soluton_vec(element_index[i - 1]) * _element ->phi_y(xi,eta,i);
            h1errorink += (u_j_y - u_h_j_y) * (u_j_y - u_h_j_y) * _element ->GaussionWeight(j);
        }
        h1error += h1errorink;
    }
    h1error /= _mesh->n_dofs();
    return sqrt(h1error);
}
#undef TEMPLATE
#else
#endif